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ABSTRACT 

^D I We present a parameter study of offset mergers between clusters of galaxies. Using the Eulerian 

^^ . hydrodynamics/A^-body code COSMOS, we simulate mergers between nonisothermal, hydrostatic clusters with a 

steep central dark matter density profile and a /3-model gas profile. We constrain global properties of the model 
clusters using observed cluster statistical relationships. We consider impact parameters between zero and five 
times the dark matter scale radius and mass ratios of 1:1 and 1:3. The morphological changes, relative velocities, 
and temperature jumps we observe agree with previous studies using the King profile for the dark matter. We 
observe a larger jump in X-ray luminosity (~ 4- lOx) than in previous work, and we argue that this increase is 
most likely a lower limit due to our spatial resolution. We emphasize that luminosity and temperature jumps due 
to mergers may have an important bearing on constraints on ft derived from the observation of hot clusters at 
^ . high redshift. Shocks are relatively weak in the cluster cores; hence they do not significantly increase the entropy 

(•~^ ' there. Instead, shocks create entropy in the outer regions, and this high-entropy gas is mixed with the core gas 

during later stages of the merger Ram pressure initiates mixing by displacing the core gas from its potential 
CSJ ■ center, causing it to become convectively unstable. The resulting convective plumes produce large-scale turbulent 

C"~~- ' motions with eddy sizes up to several 100 kpc. This turbulence is pumped by dark matter-driven oscillations in 

^D , the gravitational potential. Even after nearly a Hubble time these motions persist as subsonic turbulence in the 

cluster cores, providing 5 - 10% of the support against gravity. The dark matter oscillations are also reflected in 
the extremely long time following a merger required for the remnant to reach virial equilibrium. 

i-G ' Subject headings: galaxies: clusters: general — hydrodynamics — intergalactic medium — X-rays: galaxies 

o ' 1- introduction 

Vh ■ 
"*T* The intracluster medium (ICM) in clusters of galaxies is now understood to evolve significantly over the lifetime of a cluster due 

j^ to the complex interactions induced by mergers with other clusters. These mergers are an important part of the cluster formation 

!• process in present hierarchical models of large-scale structure. Mergers supply a substantial amount of energy to a cluster (two 

.5^ 1O'^M0 clusters colliding at ^ 1000 km s"' yield an energy ^ lO^'* erg). Thus they have been suspected of contributing to many 

K^ , different cluster phenomena, including the production of radio halos (Harris et al. 1980; Burns et al. 2000), extreme ultraviolet 

^ ' and nonthermal X-ray emission (Sarazin & Lieu 1998; Blasi & Colafrancesco 1999; Blasi 2000), the bending of narrow-angle-tail 

Ci ' and wide-angle-tail radio galaxies (Burns et al. 1994; Gomez et al. 1997; Bliton et al. 1998), and the disruption of cooling flows 

(McGlynn & Fabian 1984; Stewart et al. 1984). 

While some evidence for mergers can be obtained from optical data alone, it is the X-ray-emitting gas which provides the strongest 
evidence for mergers. The energy gained by the ICM of one cluster as it falls into the potential well of another is easily enough to 
raise its temperature to the observed values of A:r ~ 10^-10*^ K. Simulations have shown that the resulting shocks can produce 
long-lived substructure in X-ray surface brightness maps; however, they have also shown that X-ray temperature and velocity maps 
are likely to be much better diagnostics, partly because shock brightness contrasts are small when seen in projection (Roettiger, 
Burns, & Loken 1996; Ricker 1998). The Chandra and XMM satellites are producing ICM temperature maps of much higher spatial 
and spectral resolution and energy sensitivity than was previously possible (e.g., Markevitch et al. 2000; Vikhlinin, Markevitch, & 
Murray 2001; Mazzotta et al. 2001). While generally confirming that mergers produce large temperature variations in clusters, these 
new observations have raised new questions. In addition to showing high-temperature ridges characteristic of shocks, many of these 
observations show cold cluster cores moving through low density, shock-heated intracluster gas. 

The intracluster gas is repeatedly stirred during a cluster's life by merger shocks (Norman & Bryan 1998) and galaxy wakes 
(Stevens, Acreman, & Ponman 1999; Sakelliou 2000). The Reynolds number /?e associated with these motions is typically of the 
order (Sai'azin 1988) 



R^^^M^-), (1) 

where Ai is the Mach number, i is the size of the stirring agent, and A, ^ kpc is the mean free path for ions in the ICM. Merger 
shocks (A^ ^1-4) typically span lengths comparable to the size of a cluster (i ^ Mpc), so for them Re ^ 4x 10^. Galaxy motions 
(i ~ 30 kpc) are transonic relative to the ICM, yielding Re ^^ 10^. These influences on the ICM may produce turbulence; if so, it is in 
a range of Reynolds number which is accessible to present simulations. However, definitive observational detection of turbulence in 
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the ICM will require measurements of the local gas velocity field. X-ray observations are just beginning to enter the regime in which 
it will be possible to constrain these velocities using X-ray emission lines (Roettiger & Flores 2000). 

Turbulence may contribute to a number of the proposed secondary consequences of mergers, including the formation of constant- 
density gas cores, the destruction of cooling flows, and in situ acceleration of cosmic-ray particles. It can also affect the stripping 
of gas from galaxies as they orbit in the cluster potential. Unless tangled small-scale magnetic fields limit diffusion (Chandran et al. 
1999), this picture is further complicated by the breakdown of the fluid description of the ICM on kpc scales. Nevertheless, we can 
make some headway with this problem by asking whether merger-induced hydrodynamical effects on scales of tens to hundreds of 
kpc can feed energy to scales on which collisionless processes become important. It may also be possible to include collisionless 
transport effects as a correction to fluid calculations on small scales. 

Offset cluster mergers are a likely source of vorticity and turbulence in the intracluster medium. The bulk angular momentum that 
feeds this vorticity originates in tidal torques produced during the linear phase of structure formation (Hoyle 1949; Peebles 1969; 
White 1984) and should result in mergers with nonzero impact parameters. One example of such an off-center merger is Abell 754 
(Zabludoff & Zaritsky 1995; Henriksen & Markevitch 1996; Roettiger, Stone, & Mushotzky 1998; Valinia et al. 1999). 

To distinguish the effects of offset mergers from those of the many other simultaneous influences on the ICM, it is necessary to 
study them in isolation (see Schindler [2000] for a brief review). Offset mergers have been studied by only a few groups. Roettiger et 
al. (1998) and Roettiger & Flores (2000) treat them primarily in the context of attempts to model specific observed clusters. Ricker 
(1998) performed a parameter study of offset mergers between equal-mass clusters but neglected the effects of the collisionless dark 
matter component thought to dominate the mass of most clusters. Takizawa (2000) considers electron-ion nonequipartition effects 
in a parameter study of offset mergers, but these smoothed-particle hydrodynamics calculations use too few particles to adequately 
characterize merger shocks and energy transfer between the gas and dark matter. Also, published studies of controlled mergers have 
used dark matter density profiles with constant-density cores, despite the fact that simulations of hierarchical structure formation 
produce halos with cuspy central density profiles (e.g., Navarro, Frenk, & White 1997 [NFW]; Moore et al. 1998; Jing & Suto 
2000). Finally, none of these systematic studies of controlled mergers emphasizes the distribution of entropy or the development of 
turbulence. 

In this paper we describe a parameter study of offset mergers that addresses these issues. We ask the following questions. Can 
constant-density gas cores survive in collisions between clusters with steep central total density profiles? How long can the pressure 
peaks associated with these cores survive in off-center collisions? How do typical offsets affect the luminosity and temperature jumps 
expected during mergers? How is entropy generated and redistributed in mergers? How much do merger-driven motions mix the 
cluster cores and contribute to nonthermal support of the gas? 

In § 2 we describe the simulation code and isolated merger model we use in our calculations. Our single-cluster convergence study 
and the merger simulations themselves are described in § 3. In § 4, we discuss the structure of the merger remnants, the evolution of 
luminosity and temperature, and the survival of distinct pressure peaks in the simulations. Finally, § 5 presents our conclusions. We 
will discuss at greater length merger-induced turbulence and methods for detecting off-center collisions in a future paper. 

When cosmological scaling is required, in this paper we have assumed a Hubble constant Hq = lOOh km s"' Mpc"' with h = 0.6. 
Unless otherwise specified, quantities are expressed using length units of /z~' Mpc, mass units of lO'^M©, time units of Gyr, and 
temperature units of ke V. 

2. NUMERICAL METHOD 

2. 1 . Simulation Code 

For our merger simulations we have used the COSMOS A^-body/hydrodynamics code (Ricker, Dodelson, & Lamb 2000). COS- 
MOS solves the Euler equations for the intracluster medium using the piecewise-parabolic method (PPM; Colella & Woodward 
1984). The code follows the evolution of collisionless dark matter by means of the particle-mesh method (Hockney & Eastwood 
1988), solving the Poisson equation for the total gravitational potential using a full multigrid algorithm (Brant 1977). It also incor- 
porates radiative cooling of the gas, although the simulations described in this paper do not make use of this feature. All three of 
the primary code components (hydrodynamics, gravity, and dark matter) make use of static, nonuniform grids, allowing us to use 
smaller zones in regions of interest than would be possible with a uniform-grid code. Our code is fully compressible and uses a 
three-dimensional Cartesian grid, so shocks and convective effects are intrinsically included. COSMOS has been tested thoroughly 
against a suite of standard test problems, and we have already used its hydrodynamical and gravitational modules to study cluster 
mergers without dark matter (Ricker 1998). 

The combination of PPM and particle-mesh is widely used in cluster dynamics and large-scale structure simulation. PPM is to 
be contrasted with another widely-used hydrodynamical method, smoothed particle hydrodynamics (SPH; Gingold & Monaghan 
1977; Lucy 1977), a Lagrangian method which uses particles as moving interpolation centers for the hydrodynamical variables. 
We use PPM because it resolves shocks extremely well. SPH handles shocks poorly but compensates in part by achieving better 
overall spatial resolution in high-density regions. Both algorithms have been shown to give consistent results when run with the same 
cluster-formation initial conditions (Frenk et al. 1999). 

2.2. Boundary Conditions 

Clusters formed in cosmological large-scale structure calculations do not evolve in isolation. Generally they are subject to tidal 
torques in the linear phase of structure growth, and during the nonlinear phase they accrete a significant number of smaller objects, 
from galaxies to other clusters. However, we are interested here in separating out the effects of merger-generated shocks. Because 
we are studying controlled mergers in isolation, it makes sense to use outflow boundary conditions for the gas and dark matter and 
isolated boundaries for the gravitational field. As noted in § 2.3.3, this requires us to cut off the initial density profiles at a finite 
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radius, allowing particles and gas to escape through this outer boundary as the simulation progresses. Our single-cluster convergence 
study (§3.1) shows that this leads to a steepening of density profiles at large radius which is unrelated to the merger-induced changes 
in the clusters. For this reason we confine our conclusions to the region within the initial cutoff radius, focusing particularly on the 
behavior of the clusters within the innermost few core radii of their centers. 

Here we briefly describe our implementation of isolated boundary conditions; further details appear in Ricker et al. (2000). For 
hydrodynamics we use standard zero-gradient boundary conditions wherever the mass flux through the external boundary is positive. 
For edge zones with inward-directed velocities, we use zero-gradient boundary conditions for all variables except the component of 
velocity normal to the boundary; for this component we use a Dirichlet boundary condition. This procedure prevents the destabilizing 
artificial inward flows which can result when matter near the boundary begins to fall back toward the center. In the case of dark matter 
particles we discard those particles which leave the grid. To handle isolated boundaries with our Poisson solver we have developed a 
method which first obtains a zero-boundary solution, uses this to obtain boundary values for the required correction to this solution, 
then solves for the correction using the given boundary values. 

2.3. Initial Cluster Model 

2.3.1. Density Profiles 

In this paper we approach the merger problem by simulating isolated collisions between idealized model clusters. Our model clus- 
ters are initially spherically symmetric, with specified density profiles for the gas and dark matter. For the questions we address, this 
approach has two main advantages over studying clusters formed in large-scale structure simulations. We eliminate the complication 
of multiple mergers and focus on the effects of a single collision, allowing us to isolate important physical mechanisms. We also 
achieve high spatial resolution, which is critical for obtaining suitably converged results that track the length scales of interest in clus- 
ter cores. However, this approach does not consider the effects of multiple simultaneous mergers, external tidal fields, and infalling 
matter. Additionally, we neglect radiative cooling and choose our cluster parameters to ensure that the cooling timescale is larger 
than a Hubble time. We do this to avoid the computational difficulties that arise when the gas is allowed to cool very quickly (the 
so-called 'cooling catastrophe'). The effects of mergers on cooling flows present some very interesting and challenging problems, 
which we will take up in a future paper. 

In the current picture, the formation of a cluster of galaxies is a process of steady gravitational accretion and relaxation toward 
equilibrium, periodically interrupted by mergers with other partially virialized systems. These periodic interruptions destroy any 
pre-existing equilibrium. After such a merger, the cluster resumes relaxing toward a new equilibrium state, continuing until it is 
again interrupted. Previous calculations have shown that this relaxation time freiax is long compared to the likely merger interval fmerge 
(Roettiger et al. 1993; Schindler & Miiller 1993; Pearce, Thomas, & Couchman 1994; Ricker 1998), meaning that mergers usually 
involve clusters which are not fully virialized. However, studying the case of large fmerge/freiax is a necessary preliminary to study of 
the case in which a merger perturbs a partially virialized initial state. In this paper we are interested in understanding the physical 
mechanisms important in individual mergers rather than simulating a specific observed cluster 

Accordingly, we begin with an initial model that represents our best current knowledge of the structure of fully virialized, nonin- 
teracting clusters. We assume functional forms for the gas density and total density profiles, then set the dark matter profile, the gas 
temperature profile, and the dark matter velocity dispersion profile in such a way as to establish hydrostatic equilibrium. The gas 
density profile is taken from the /3-model (Cavaliere & Fusco-Femiano 1976), 

-3/3/2 
/ r \ 
1 + 



Ps(r) = PgQ 



(2) 



where the central density pgo, core radius r^, and asymptotic slope parameter f3 are fitting parameters. When used with the assumption 
of isothermality, the /3-model with f3 w 2/3 produces an X-ray surface brightness profile which is consistent with those observed in 
many clusters without cooling flows (Jones & Forman 1984). In this paper we take f3 = 2/3 as a representative value. The /3-model 
is similar in shape to the nonsingular isothermal sphere profile (Binney & Tremaine 1987), which is a solution of the gas hydrostatic 
equation under the assumption of isothermality. 

We take guidance for the form of the dark matter density profile from models of hierarchical large-scale structure formation. Using 
A^-body simulations, NEW have shown that, in such models, dark matter halos on scales from individual galaxies to rich clusters 
form with density profiles given by 

21-1 



P(r) = Ps 



r 

1 + - 
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(3) 



where p, and r, are, respectively, a scaling density and radius, both of which depend upon the mass of the halo in question. The mass 
dependency in these simulations arises because of the differing redshift of formation of halos with differing masses, with the sense 
of the dependency being that less massive halos, which form earlier, have higher central densities and are thus more concentrated. 
The exact shape of this universal density profile has been called into question by Moore et al. ( 1 998), who favor a steeper profile near 
the centers of halos, arguing that previous simulations have insufficiently resolved the halos. Nevertheless, the clear prediction of 
these models is that self-similar, purely gravitational evolution leads to a steeper density profile than that predicted by the /3-model. 
The addition of collisional, self-gravitating gas most likely modifies this profile somewhat, though if dark matter dominates as we 
expect it to, the total density profile should still be steeper than that of the gas near the center. Therefore, we choose a total density 
profile using the NEW functional form. The dark matter density profile is given by the difference between the total profile and the 
gas profile: 

Pdm(r) = p(r)-pg(r) . (4) 
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This can only be an approximation to the true situation, because the NFW profile has been established by A^-body simulations which 
do not include gasdynamics. If the NFW profile were to be maintained for the total density upon including gasdynamics, the shape of 
the dark matter profile would need to be different enough to exactly compensate for the /3-model shape of the gas. However, because 
the dark matter dominates the total density, particularly in the center of the cluster where the NFW profile continues to rise while the 
/3-model tends to a constant value, this is most likely an acceptable approximation in the inner regions of non-cooling-flow clusters. 

2.3.2. Temperature and Velocity Dispersion Profiles 

We choose the gas temperature and dark matter velocity dispersion profiles from the equations of hydrostatic equilibrium for the 
two matter components. A more complete discussion will be given in a subsequent paper (Ricker & Sarazin 2001, in preparation). 
With the total density profile given by equation (3), the total mass M enclosed within radius r is given by 

X 



M(r) = 47rp.,rj 



ln(l+x)- 



l+x 



where x = rjr^. The hydrostatic equation for the gas temperature profile r(r) is 



d_ 
dr 



Pg{r)kT(r) 
^niH 



GM(r) 



Pg(f) , 



where ^iniH is the average mass per particle in the gas. It is useful to define a characteristic temperature scale 7} as 

■* 7 



Then the dimensionless temperature function 7'(x, T,) is defined by 



r(xr,) 



If we define rj as the ratio of the total mass scale length to the gas core radius, 

r. 
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then the gas density can be written as 

where we have used our assumption that (3 = 2/3 

equation becomes 

d 
dx 



p,(r) = pjo(l + r/V)-\ (10) 

Using equations (5) through (10), the dimensionless form of the hydrostatic 
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(11) 



In the case of /3 = 2/3 which we are considering here, the solution to this equation (with the boundary conditions discussed below) 
can be written in closed form in terms of Clausen's integral. However, for the purpose of initializing our simulations it is more 
straightforward to solve this equation numerically. 

The outer boundary condition for equation (11) depends on the pressure at large radii in the cluster One class of solutions would 
have the pressure approach a constant value at large radii to maintain pressure equilibrium with surrounding intercluster gas. Another 
might have the pressure go to zero at a finite radius, which would represent the outermost extent of the hot gas in the cluster This 
condition might approximate a cluster with a strong accretion shock at its outer radius, with the gas pressure being very small outside 
of this radius. However, the overdensity in clusters is large (> 10^ in the core), so the outer pressure is likely to be very much smaller 
than the central pressure. Thus, we adopt the boundary conditions that apply to an isolated cluster that extends to large radii, so that 
the pressure approaches zero at large radii. Assuming the adopted form for the gas density profile (eq. [2]), this implies that f{x) -^ 
as X ^ oo. With this boundary condition, the temperature profile requires no additional scaling parameters beyond those supplied by 
the assumed density profiles, and it is independent of the normalization of the gas profile (p^o)- 

We constrain the allowable values of j] by noting that convective stability requires 

??>0.71, (12) 

while in order for the central cooling time to be larger than a Hubble time (to avoid producing a cooling flow), we must have 

-3 \ P" \ I r, \ ( pgo 



V ^ 9.6/;" 



(13) 



lOi5M0(/!-i Mpc)-V V/z-i Mpcy \W^MQ(h-^ Mpc)- 
(Ricker & Sarazin 2001). By experimenting with typical values of p,, pgo, and r,, we have found 77 = 2 to be a safe value given the 
resolutions we obtain in the simulations described in § 3. This is comparable to the sizes of the numerical core radii seen in some 
large-scale structure calculations (e.g., NFW). 

The resulting gas temperature profile (scaled to the central temperature) is shown as the dotted line in Figure 2e (below). For 
X < 2/77 the temperature drops slowly (roughly as x~°°^), then flattens out and eventually tends toward an asymptotic ln(x)/x behavior. 
By the radius x = 20/77, the temperature has dropped to 25% of its central value. 

We follow a similar procedure to determine the dark matter velocity dispersion profile. We assume the dark matter velocity 
components to be isotropic and Gaussian distributed with velocity dispersion <7^(r). The virial equation for the dark matter particles 
is then 

-J- [pAm(r)a\r)\ = ^pdm(r) . (14) 

dr r^ 
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For the dark matter the density profile we use contains p,, pgo, Vs, and r^; hence no simple, natural scaling for a is available (as in 
eq. [7] below). Therefore we choose the velocity dispersion at the cutoff radius r = R,a\ = cP'ifi), and constrain its value using the 
assumption of virial equilibrium (eq. [26]). As with the gas hydrostatic equation, we numerically integrate the resulting dimensionless 
virial equation. The result {for R/rs =10, which we motivate in the next section) appears as a dotted line in Figure 2b. The velocity 
dispersion rises from its central value to about twice this value at x w 1, then drops slowly back to its central value at r = /?. The 
difference in behavior at small radii in comparison to the gas temperature is due to the differing behaviors of the dark matter and gas 
densities near the center. The dark matter requires a 'colder' center to maintain its cuspy density profile, while the gas requires a 
hotter center in order to maintain its constant-density core. 

2.3.3. Scaling the Cluster Model 

By choosing the dark matter and gas density profiles, we introduce as model parameters the density scales ps and pgo, and the 
length scales r, and re- Because both density laws (for typical values of j3) yield divergent total masses for r ^ oo, we introduce a 
cutoff in the density profile at r = /?, making R an additional parameter. We examine the effects of this cutoff in the single-cluster test 
calculations discussed in § 3.1. The velocity dispersion profile introduces the additional 'thermal' parameter ctr, yielding a total of 
six parameters. (The scaling parameter T, for the gas temperature profile is determined by p, and r, via equation [7].) 

While we could examine the evolution of clusters with chosen values of the above parameters, it is more useful to select values 
for the total mass M, X-ray luminosity Lx, emission-weighted temperature Tx, and gas fraction Fg, expressing the density, thermal, 
and length scales in terms of these quantities. This enables us more readily to compare our results with real clusters. In addition, it 
permits us to use observed correlations between these global quantities to reduce the dimensionality of the parameter space we must 
study. 

We begin by defining the dimensionless density and dark matter velocity dispersion profiles: 

p{xr,;) 
p(x) = (15) 

pjxr,) 

Pg(x,fg,7^)=^^S^^ (16) 

Ps 

(7(xrs) 
a(x,aR)= -. (17) 

The dimensionless gas temperature profile is defined in equation (8). Here we define 

/, ^ ^ (18) 

Ps 

R 

c=-. (19) 

If we take R to be the virial radius r20(), c becomes the familiar halo concentration parameter. While NEW find in their simulations 
that c has a weak dependence on halo mass, we adopt the simplification of a constant value c = 10, consistent with their results for 
cluster-mass halos. Together with our constraint on rj (§ 2.3.2), this reduces the number of free parameters to four. 
We can now write expressions for our desired quantities M, Lx, Tx, and Fg: 

M = ATTp,rl / dxx^pix) (20) 

Jo 

2.3 r ^^JIt, 



Lx=A^pirl dxx'pg(xJg,riYK[T(x)] (21) 

Jo 

Tx = ^-^p^ f dxx^pg{xJg,rifK[T{xW{x) (22) 

l^x Jo 

L dxx^pg(x,fg,Ti) 
Fg = -"^ '^: '■'^' " . (23) 

Jg dxx^pix) 

Here A(r) is the plasma emissivity function; we use the emissivity function of the MEKAL model in XSPEC (Mewe, Kaastra, & 
Liedahl 1995) with half-solar abundances. With the chosen values of ry (§ 2.3.2) and c, only three of these equations are independent; 
one of the four quantities on the left-hand side, say Tx, depends on the others. Therefore, given M, Lx, Fg, and our constraints on t] 
and c, we can numerically solve the three independent equations for p,, r,, and fg. We also have cts as a remaining free parameter. 

We can further reduce the number of free parameters by considering 'typical' virialized clusters. Because M and Fg are obser- 
vationally the most poorly determined quantities, it makes sense to allow them to vary. We constrain Lx in terms of M using the 
observed luminosity-temperature relation and the virial mass-temperature relation: 

Lx = KlT^' (24) 

M = KmT^" . (25) 

The observed exponent of the Lx- Tx relation lies in the range ^^2.5 to ^-^3.0 (Markevitch 1998; Arnaud & Evrard 1999; Allen 
& Fabian 1999) when corrected for effects due to cooling flows. We have used the bolometric relation determined by Markevitch 
(1998) from a sample of 35 nearby clusters. In our notation Markevitch's relation corresponds to Ki = 7.62 x lO'*^ erg s"' keV"^^'* 
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and fl£ = 2.64. Mushotzky and Scharf (1997) show the cluster Lx - Tx relation to be nearly independent of redshift up to z ^ 0.4; 
hence this is likely to be a good approximation even for clusters at moderate redshifts. 

The virial mass-temperature relation we use is taken from Evrard, Metzler, & Navarro (1996), who find in an ensemble of 58 
simulated clusters that equation (25) with Km = 7.02 x IO'-'Mq keV"-'''^ and gm = 3/2 predicts cluster masses to within a standard 
deviation of 15%, even for fairly irregular clusters. 

Finally, we constrain (Jr, the scaling factor for the dark matter velocity dispersion, using the assumption of virial equilibrium 
within r = R. We also assume that dark matter particle orbits at r = R are nearly isotropic. This yields the constraint 

, GM 
4 = ^. (26) 

In setting up the initial particle distribution for each cluster, we choose the x-, y-, and z-components of the velocity for each particle 
from a Gaussian distribution with dispersion <J^(r). 

2.4. Choice of Parameters 
2.4.1. Model Clusters 

With the assumptions discussed above, the properties of each cluster are completely determined by the assumed mass M and gas 
fraction Fg. We chose gas fractions which fell within the range of values permitted by X-ray cluster observations (David, Jones, & 
Forman 1995; White & Fabian 1995) and which yielded central cooling times fcooi.o larger than a Hubble time, //^ '. 

Table 1 presents the parameter values for the two cluster models studied in the collision runs. Any other cluster properties can be 
derived from these values using the scaling relations discussed above. The temperatures of the two clusters were chosen to fall in the 
range of 'typical' rich clusters seen at low redshift (e.g.. Edge et al. 1990; Edge, Stewart, & Fabian 1992). We then chose masses 
using the virial mass-temperature relation as described in section 2.3.3. The values of the scaling radius r, and central density (given 
in Table 1 in terms of the central electron number density, neo) were constrained as described in § 2.3.3. Table 1 also shows the 
calculated values of the sound-crossing time fsc = R/csO and the 2-10 keV X-ray luminosity L^~'". The two clusters, which we label 
'A' and 'B', have a mass ratio of approximately 2.8. Among real nearby clusters, our simulated clusters A and B are most nearly 
similar to Virgo and AWM7, respectively, although both of these real clusters have weak cooling flows. 

2.4.2. Initial Merger Kinematics 

Our collision model, depicted in Figure 1, is appropriate for two clusters that have fallen together from a large distance with 
nonzero angular momentum. This model introduces three new parameters: the initial separation d, the initial relative velocity 
V, and the impact parameter b. Because clusters of given masses do not all merge at the same time or with the same total angular 
momentum, we have some freedom in choosing the values of these parameters. However, only certain ranges of values are physically 
well-motivated. Here we discuss the scaling of our choices with the properties of our model clusters. 

For computational reasons, it makes sense to fix the separation d at the start of the calculation, balancing the need to save compu- 
tational time in following the relatively uninteresting pre-merger cluster evolution with the need to start the clusters from a virialized 
state. In each case we set 



d=^(Ri+R2r + b\ (27) 

where Ri and R2 are the cutoff radii of the two clusters. We assume the clusters to have fallen to this separation from a greater 
separation d{) w l.6d, at which their relative radial velocity is zero. The value of do corresponds to the turn-around distance for the 
clusters (Sarazin 2001), assuming that the cluster merger (the first core crossing) occurs when the age of the Universe is ^1 1 Gyr 
(the present age for il = 1 and h = 0.6). The exact value of do does not affect the total energetics much as long as the relative velocity 
of the clusters at closest approach is close to the value expected for free-fall from infinity. 

For the purpose of computing the initial relative velocity, we approximate the two clusters as point masses. Comparing the 
angular momentum and energy predicted with this approximation to the actual value determined from the simulations shows that 
this estimate is good to within 1% because of the highly concentrated nature of the cluster density profiles. We also neglect the 
clusters' spin angular momenta. At the separation do, the clusters are assumed to have zero relative radial velocity; hence their 
angular momentum and energy are 

Jq « mvodo (28) 

1 2 GM1M2 
^°"2'"^« — dT' 

Table 1 
Scaling Parameters for Model Clusters 



Cluster 


M 

(IO'^Mq) 


Tx 
(keV) 


F, 


(/j-i kpc) 


rieO 

(10-^ cm-3) 


(Gyr) 


fcool.O 

(Gyr) 


Lx(2-10keV) 
(lO'^^ergs-i) 




A 
B 


0.199 
0.562 


2.0 
4.0 


0.125 
0.174 


118 
169 


3.95 
5.31 


1.89 
1.93 


17.7 
19.7 


0.117 
1.25 
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Fig. 1 . — Schematic of the initial conditions used for cluster merger calculations. The clusters have masses Mi and M2 and cutoff radii Ri and R2- The initial 
separation of the cluster centers is d, and the impact parameter (j'-direction) is h; the collision plane is the x>'-plane. The initial relative velocity is v; the simulation 
takes place in the center-of-mass frame. 



where their reduced mass is 



M1M2 

M1+M2 



(29) 



and vo is their relative azimuthal velocity. At the separation d, the relative velocity v is perpendicular to the direction of b, so we can 
write 



J « mvb 

1 2 
E w —mv 

2 



(30) 



GM1M2 



Conserving angular momentum and energy, we eliminate vo and find 



1 1 



v^ = 2G(Mi+M2) 

d do 



do 



H-l 



(31) 



We estimate the range of impact parameters b to study by using the linear-theory result for the dimensionless spin of dark-matter 
halos. The spin parameter A is defined as (Peebles 1969) 

Here J is the angular momentum of the halo, E is its total energy, and M is its mass. Recently, Sugerman, Summers, and 
Kamionkowski (2000) have performed a detailed comparison of linear-theory predictions to actual angular momenta of galaxies 
formed in cosmological A^-body/hydro calculations. These simulations did not include cooling or star formation, so at the upper end 
of the mass range they studied their results should carry over to the clusters we are simulating. They find, in agreement with White 
(1984), that linear theory overpredicts the final angular momentum of galaxies by roughly a factor of three, with a large (~ 50%) 
dispersion in the ratio of the linear-theory prediction to the actual value. The importance of nonlinear effects on the spin parameter 
is less clear, as uncertainties in their group-finding algorithm strongly affect the determination of J, E, and M for each galaxy and 
produce a large scatter in A. However, given this uncertainty, the mean value of A they obtain is about 0.05, close to the linear-theory 
result. 

We take the halo to be the final merger remnant. Its final total angular momentum is the sum of the spin angular momenta of the 
two subclusters plus the orbital angular momentum /orb- We estimate the initial spins of the subcluster by applying equation (32) to 
each of them, and we assume that the initial spins are correlated. The orbital angular momentum is then the difference between the 
final angular momentum of the merger remnant, and the total spin angular momenta of the subclusters (Sarazin 2001). Substituting 
equations (28)-(31) into equation (32) and solving for the impact parameter, we find that 

b^M'-^y'(l-4-Y''muM2). (33) 



l-jj f(MuM2). 



Here, the function /(Mi ,M2) corrects for the internal angular momenta and energy of the subclusters, and can be written as 



f(MuM2) = 



(Mi+M2r 



(Ml +M2f'^ 



3/2 



(34) 



but it only depends on the ratio (M^/My) of the smaller to larger mass of the two subclusters. It varies between 4(2^/^- 1)^/^ « 
1 .80 < /(Ml ,M2) < (5/3)^/^ « 2. 15, so that /(Mi ,M2) « 2. In our merger models, d « 20rj^max ('"^.max is the larger of the two scale 
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radii of the clusters), andJo w l.6d. If we assume A «0.05 and approximate /(Mi, M2) ~2, we find that the typical impact parameter 
is fe w Srj.niax- Since this argument is statistical and a range of values are expected in mergers, we consider values of fe = (0, 2, 5) rj.niax 
in our simulations. 

Table 2 presents the parameter values used in the collision runs. The simulations are labeled 'CI' to 'C6'; runs C1-C3 are equal- 
mass collisions using cluster A, while runs C4-C6 use clusters A and B. Each run is followed for approximately 15 Gyr, more than 
the age of the Universe for 51 = 1 and h = 0.6. For the hydrodynamics and gravitational potential calculations, we use a nonuniform 
mesh with 256 x 128^ zones covering a volume measuring approximately 6/?n,ax x (A-Rmax) , where Rmax is the larger of Ri and R2. 
The minimum zone spacing is approximately equal to one-seventh of the larger of the two dark matter scale radii, or ^ 17/;"' kpc 
for runs C1-C3 and ^ 25/z~' kpc for runs C4-C6. For the dark matter, each run also uses 128^ particles of equal mass, distributed 
between the merging clusters in proportion to their total dark mass. The results of these runs are described in §§ 3.2 and 3.3. In § 3.1 
we first discuss the results of a convergence study involving a single cluster, which we used to set the grid and particle parameters 
used in the collision runs. 

3. NUMERICAL SIMULATIONS 
3.1. Single-Cluster Convergence Study 

We carried out several calculations of the motion of a single cluster with different spatial resolutions and numbers of particle in 
order to evaluate numerical effects due to our boundary conditions and finite-radius cutoff. These runs also enabled us to choose 
grid and particle parameters for the collision calculations. These runs used the cluster parameters appropriate to cluster A (see § 2.4) 
together with the numerical parameters listed in Table 3. Six runs were performed with two different grid spacings and three different 
numbers of particles. One additional run was performed at a grid spacing and particle number slightly better than the collision runs 
to be described later. The computational volume measured 7.08 x 4.72 x 4.72(/z"' Mpc)-', and the cluster started at the position 
(2.36,2.36,2.36) /z~^Mpc with an initial velocity of 0.436 /i~'Mpc Gyr~' in the x-direction. In each run the cluster was permitted to 
move for 5 Gyr, yielding an expected final displacement of 2. 1 8 /;"' Mpc. 

Figure 2 compares final average profiles of various quantities in the single-cluster tests to the initial profiles. Included are the dark 
matter density and velocity dispersion, gas density, enclosed gas fraction, gas temperature and specific entropy, and the gas radial 
Mach number. We also plot the ratio of the gas pressure gradient to the gravitational force per unit volume on the gas, which shows 
the extent to which the gas is in hydrostatic equilibrium. The results show the effects of finite zone spacing and particle number on 
small scales and outflow boundary conditions on large scales. For gridded (gas) quantities, each profile was generated by averaging 
profiles interpolated along 1000 different directions from the gas centroid. The radial sample spacing was taken to be equal to the 
size of the smallest zones in the calculation. For particle (dark matter) quantities, particles were binned in radial bins about the dark 
matter centroid; then the particle count or the average value of the velocity dispersion in each bin was used. The bin spacing used 
was the same as the sample spacing used for the gas profiles. With the exception of the gas fraction, the pressure-potential ratio, and 
the radial Mach number, all quantities are scaled in terms of the initial values of the appropriate scaling parameters, such as the core 
radius or central gas density. 

The dark matter density profile (Figure 2a) shows little sensitivity to the number of particles or zone spacing for the ranges of 
these parameters considered here. Resolution effects appear in the innermost part of the cluster; the density in the innermost radial 
bin tends to be some 15% below the initial value. This demonstrates the tendency of particle-mesh force smoothing to produce an 
artificial constant-density core on scales smaller than twice the mesh spacing. Thus in order to observe the turnover to r~' behavior 
at small scales characteristic of the NFW profile, the mesh spacing should be no larger than about 1/4 the scale radius r,. For r > Ir,, 
the escape of particles through the outer boundary of the cluster makes the resulting final profile slightly steeper than its initial r~^ 

Table 2 
Cluster Parameters Used for Collision Runs 





Cluster 


Mass 








V 


4 4. 


t 




Identifier 


Models 


Ratio 


b/ 


'.v.max 


(h- 


' Mpc Gyr-i) 


(h-^ Mpc)'' 


A ■ /r • ^''^ 




CI 


AA 


1 : 1 









0.436 


7.10 4.73 


4.73 


0.144 


C2 


AA 


1 : 1 




2 




0.436 


7.10 4.97 


4.73 


0.144 


C3 


AA 


1 : 1 




5 




0.436 


7.10 5.32 


4.73 


0.144 


C4 


AB 


1:3 









0.548 


8.63 6.78 


6.78 


0.213 


C5 


AB 


1:3 




2 




0.548 


8.63 7.12 


6.78 


0.213 


C6 


AB 


1:3 




5 




0.548 


8.63 7.63 


6.78 


0.213 



^^s.max and rv.min ^6 the larger and smaller of the two scale radii, respectively. 
^£x, £y, and f ^ are the physical dimensions of the computational grid. 
'^ Amin is the smallest zone spacing in any direction. 
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shape. Runs SI -S3 display a somewhat shallower profile between 2rs and Sr, than that in runs S4-S7. 

The dark matter velocity dispersion profile (Figure 2b) also reflects this evaporative cooling: within r w r, in run S7, cP' is about 
9% lower than its initial value, with roughly the same r dependence. Outside r w r, the final profile drops off more steeply, whereas 
outside r = R ii appears to rise again due to the radial motion of escaping particles. (We estimate cr^(r) by computing (v?) - (v,)^ 
for each velocity component /, averaging within radial shells, then sum the results for x, y, and z and divide by 3. This procedure 
assumes isotropic particle velocity dispersions and does not subtract out radial expansion or contraction.) The central dispersion in 
runs SI -S3 is substantially lower than in run S7, but runs S4-S6 agree with S7. 

The lack of an external confining pressure produces expansion of the gas at the outer boundary of the cluster. This has a major 
effect on the gas density profile (Figure 2c), causing it to become steeper at large radii. The profile also becomes slightly shallower 
at small radii. However, unlike the dark matter profile, once the gas core radius is resolved (runs S4-S7), the gas profile converges to 
a central density slightly below the initial value. Between r ^ Vc and r ~ Avc it agrees with the initial profile, but the effective core 
radius (defined as that radius at which the density drops to one-half its central value) at the end of the simulation is ^ 2% smaller 
than the initial value, so the overall shape at small radii differs slightly from a /3-model. When the core is not resolved (runs SI -S3), 
the behavior is more like the dark matter density, with the innermost density point well below its initial value. Again, the number of 
particles makes little difference for the range of values considered here. 

The different behaviors of the gas and dark matter are apparent from the enclosed gas fraction profile (Figure 2d). The enclosed 
gas fraction /^gas(< r) is computed from the average gas and dark matter densities via the relation 

Jq r-^drptotin 

In all cases the gas fraction at a given radius is smaller than its initial value for radii larger than about two zones. Within two zones of 
the center the gas fraction flattens out to a constant value rather than dropping to zero because of the artificial core in the dark matter 
profile. The final gas fraction rises more slowly than its initial profile, reaching the initial total value of 12.5% only at r ^ 2R. 

The spreading of the gas due to the lack of external pressure confinement leads to a steepening of the temperature profile (Figure 
2e) and the advection of entropy out of the outer parts of the cluster (Figure 2f). Whereas the initial temperature profile drops to 
~ 0.4 times its central value at r = R, by the end of the single-cluster calculation it has dropped to about 0.2 times its central value. 
In runs SI -S3 the gas core is not resolved, and the innermost temperature point is about two-thirds its initial value at that radius. For 
runs S4-S7 the final temperature agrees much better with its initial profile, although run S7 does not reproduce the upturn within 
^ O.Src present in the initial conditions. The two sets of core-resolved runs also disagree slightly between r ^ Vc and r ^ 3rc, with 
runs S4-S6 agreeing better with the initial profile. All of the simulated profiles agree with each other beyond about five core radii. 
Similar results can be seen in the specific entropy profile. Note that the specific entropy is measured relative to the central value 
given the imposed density and temperature profiles. The resulting imposed entropy profile passes through a minimum near r = O.Sr^, 
so at smaller radii the gas would be convectively unstable if our calculations were to use smaller zones. 

The radial Mach number (Figure 2g) and pressure support profiles (Figure 2h) show that the innermost parts of the cluster approx- 
imate hydrostatic equilibrium much better as grid resolution improves. Both quantities are formed from angle-averaged values; thus 
the radial Mach number is computed using the average radial velocity and the average temperature. When the core is not resolved, 
the innermost radial velocity point is as much as 0.2 times the sound speed, and the pressure provides about 10% of the gravitational 
force there. In runs S4-S7 the maximum Mach number inside r = lOr^ is less than 5%, and at the innermost point the pressure 
contributes 30% of the force due to gravity. The pressure support ratio increases immediately to a value slightly larger than unity 
outside the center, indicating that the outer parts of the cluster are nearly in hydrostatic equilibrium but are still slowly expanding. 
In all cases the gas outside r = lOr^ is in nearly free expansion, with the radial velocity increasing to a sonic point at r ^ 25re. The 
material inside this sonic point is slowly decelerating and will eventually fall back onto the cluster. 

In summary, the single-cluster test calculations show that the values of particle number and zone spacing used in the collision runs 
are adequate to produce a control cluster which remains in nearly hydrostatic equilibrium. The code is especially well-converged 

Table 3 
Mesh/Particle Parameters Used for Single-Cluster Runs 



Identifier 


Mesh Size 


No. of Particles 


A„,i„/r/ 




SI 


64 X 32 X 32 


65536 


0.80 


S2 


64 X 32 X 32 


524288 


0.80 


S3 


64 X 32 X 32 


1048576 


0.80 


S4 


128 X 64 X 64 


65536 


0.25 


S5 


128 X 64 X 64 


524288 


0.25 


S6 


128 X 64 X 64 


1048576 


0.25 


S7 


256 X 128 X 128 


1048576 


0.12 



is the smallest zone spacing in any direction. 
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Fig. 2. — Average profiles of dark matter and gas quantities at ? = 5 Gyr in the single-cluster test runs. Parameters used in the test runs are given in Table 3. In all 
cases except (d), (g), and (h), quantities are scaled to their initial characteristic values, (a) Dark matter density, (b) Dark matter velocity dispersion, (c) Gas density, 
(d) Enclosed gas fraction, (e) Gas temperature, (f) Gas specific entropy, (g) Gas radial Mach number, (h) Ratio of pressure force to gravitational force on the gas. 

with respect to particle number for the value (~ 10^) used in the collision runs, and the quantities of interest have converged down to 
radii of about two zones for the zone spacing used (~ rc/3 .5). The accuracy of the converged solution is affected by the particle-mesh 
smoothing length on small scales and by the lack of an external confining pressure on large scales. These effects lead to an artificial 
core radius of 1-2 zones in the dark matter and a nearly free expansion of the gas and dark matter outside the initial density cutoff 
which steepens the density and temperature profiles near the cutoff. 



3.2. Equal-Mass Mergers 

In this section we describe the simulations of mergers between clusters of equal mass (runs C1-C3). Figure 3 shows as functions 
of time the projection of the gas and dark matter densities along the z-axis, which is perpendicular to the plane of the collision in the 
offset runs. We refer to the projected densities as Sg and Sdm, respectively. Figure 4 shows the 2-10 keV X-ray surface brightness 
Sx and X-ray emission-weighted temperature Tx, also projected along the z-axis. Figure 5 shows the total luminosity Lx and average 
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emission-weighted temperature for the entire computational grid as functions of time. The zero of time in each plot has been chosen 
to be the instant of maximum luminosity (in each case approximately 3.5 Gyr or 1 .8fsc after the start of the calculation). The abscissa 
of the luminosity and temperature plots is also scaled to the sound-crossing time fjc- 

As the outer regions of the two clusters begin to interact early in the merger, a pair of roughly planar shocks forms at an angle to 
the collision axis, moving in opposite directions. By 0.5 Gyr prior to maximum luminosity, these planar shocks have passed through 
the cluster cores and are on their way out of the interaction region. They are readily apparent in the projected gas mass plot (Figure 3) 
and are barely visible in the surface brightness contours (Figure 4). Between the clusters the temperature is higher than the pre-shock 
value by about a factor of three, and the material in this region is expelled from the collision axis at about 1/2 the initial relative 
velocity. The material between the clusters is adiabatically compressed after the passage of the shocks. In previous gas-only merger 
calculations (Ricker 1998), the planar shocks later developed into spiral shocks which acted to dissipate the angular momentum of 
the cluster cores into large-scale rotational motion. With the addition of collisionless dark matter, however, the planar shocks play 
less of a role in the later stages of the collision, and we see more small-scale rotation as early as 3^ Gyr after maximum brightness 
and continuing for at least 3 Gyr beyond that. 

Although we have not attempted to reproduce any specific clusters with these calculations, we note that the simulations do resemble 
some merging clusters with particularly simple geometries. The temperatures and X-ray surface brightnesses in the early phases of the 
zero-offset, equal mass mergers strongly resemble the observed ROSAT image and ASCA temperature map of the cluster surrounding 
Cygnus-A (Markevitch et al. 1999). The offset 1:1 cases, seen just after first core passage, resemble the structure of Abell 3395 
(Markevitch et al. 1998), which appears to be observed from above the merger plane. 

The merging clusters undergo a dramatic increase in luminosity and temperature as the cluster cores collide. The luminosity jump 
ranges from a factor of 4 for the most off-center case to a factor of 10 for the head-on case. The average temperature increases by a 
factor of two to three. In the projected mass and X-ray maps, the period of core collision corresponds to the formation of a strong, 
crudely ellipsoidal shock which dissipates some of the initial kinetic energy of the cores into thermal energy in the surrounding gas. 
The greater increase in total luminosity and average temperature for the head-on case results from the fact that the dark matter cores, 
whose gas content produces the majority of the X-ray emission, pass directly through one another and the adiabatically heated gas 
between them. In the offset cases this hot gas is driven ahead of the cluster cores as they slip past one another, producing an S-shaped 
region of high temperature. In all cases, however, the dark matter cores continue to have associated peaks in the gas density even 
after they pass one another, a significant departure from the behavior in the gas-only collisions considered by Ricker (1998). Here 
the central potential is dominated by collisionless dark matter, which dissipates energy in a relatively inefficient manner compared 
to the gas during core passage. Thus the dark matter cores, which are able to pass through one another even in the head-on case, 
drag gas density peaks along with them. The extra oscillation due to the dark matter allows the dark matter to 'pump' the gas, 
keeping the system disrupted longer than would be expected if both matter components were collisional. Offsets on the order of a 
few hundred kpc are sometimes observed between the dark matter and gas density maxima, particularly as the cluster cores reach 
maximum separation after their first interaction. 

The first period of enhanced luminosity and temperature lasts for approximately one sound-crossing time. After their close 
interaction, the cores recede from one another, causing the temperature to drop to about 80% of its initial value. By the time 
the cores reach their maximum separation, the luminosity decreases much more, reaching about 15% of its initial value in the most 
offset run and about 30% in the head-on run. As one might expect due to the relative efficiency of gas dissipation in the three cases. 



b=Q 



b=2r. 



b-5r. 




-D.5 G^r 
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-1-2 Gyr 



-1-3.5 Gyr 



Fig. 3. — Surface density maps for dark matter (colormap) and gas (contours) in tlie 1:1 collision iijns (C1-C3), projected along the z-axis. A small region 
[~ (A"' Mpc)^] of the projected simulation volume is shown. Logarithmic color scaling for the dark matter surface density Sjn, is indicated by the key to the right 
of the figure; units are 10'^ Mq (ft^'Mpc)"^. Contours of gas surface density Sgas are spaced by a factor of two, with the outermost contour having the value 10"^ 
in the same units. The fiducial bar in the lower left comer is 0.5/7"' Mpc long. For each run, t = corresponds to the time at which the luminosity of the system 
reaches its maximum value. 
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Fig. 4. — 2-10 keV emission-weighted temperature 73c (colormap) and X-ray surface briglitness Sx (contours) in tlie 1:1 collision runs (C1-C3), as seen along 
the z-axis. A small region (~ (ft"' Mpc)-) of the projected simulation volume is shown. Linear color scahng for Tx is indicated by the key to the right of the figure; 
units are keV. Surface brightness contours are spaced by a factor of three, with the outermost contour having the value 6.3 X 10" erg s" cm"^. The fiducial bar in 
the lower left comer is 0.5/t~' Mpc long. For each run, t = conesponds to the time at which the luminosity of the system reaches its maximum value. 




Time (i-f^J/i,, 

Fig. 5. — 2-10 keV luminosity Lx and average emission-weighted temperature Tx as functions of time t in the 1:1 collision runs (C1-C3). Luminosity and 
temperature are each scaled to their initial values. Time is offset from the time of peak luminosity fn,ax and scaled to the sound-crossing time fsc ■ 
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the clusters require longer to reach maximum separation in run C3; this occurs about fsc after core passage in this run, compared to 
only 1/2 sound-crossing time for run CI. In each case, a second core passage occurs between 1.5 and 2.5 sound-crossing times after 
the first, producing smaller brightness and temperature jumps. Following this second passage, the total luminosity drops to between 
30 and 70 percent of its initial value, while the average temperature drops to 1.4 times its initial value. The cores continue to oscillate, 
producing minor variations in the luminosity and temperature, but the values of these quantities change very little for the remainder 
of the calculation. During this time several weaker sets of shocks are generated by core passages, as is apparent at 3.5 Gyr after the 
first core passage in the projected mass and X-ray plots (Figures 3 and 4). We discuss the evolution of luminosity and temperature in 
more detail in § 4.2. 

After the cores have coalesced, angular momentum continues to be redistributed in the gas by turbulent eddies which are fed by 
material falling back in the wake of the ellipsoidal shock. The remnant requires several billion years to reach equilibrium, and indeed 
has not quite reached equilibrium by f = 15 Gyr, the end of the calculation (and more than a Hubble time for h = 0.6). We consider 
the properties of the merger remnant at / = 15 Gyr in § 4. 1 . 

3.3. Mergers with a 1:3 Mass Ratio 

In this section we describe the simulations of mergers between clusters with a mass ratio of approximately 1:3 (runs C4-C6). 
Figure 6 shows mass maps projected along the z-axis (Eg and Edm) for these runs at several different times. Figure 7 shows the 
2-10 keV X-ray surface brightness Sx and X-ray emission-weighted temperature Tx, also projected along the z-axis. Figure 8 shows 
the total luminosity Lx and average emission-weighted temperature as functions of time. In these figures the more massive cluster 
initially approaches from the right. As for the 1 : 1 merger runs, the zero of time in each plot has been chosen to be the instant 
of maximum luminosity. The abscissa of the luminosity and temperature plots is also scaled to the sound-crossing time fsc of the 
more massive cluster. Numerically, the values of the time of maximum luminosity and the sound-crossing time are similar to those 
obtained in the 1:1 collisions (^ 3.5 Gyr and 1.9 Gyr, respectively). 

In the 1:3 runs the initial pair of shocks is not planar; instead, the shock sweeping through the less massive cluster curves around 
its core, and the adiabatically compressed region between the clusters takes on an arc shape. For nonzero impact parameters, the 
shock being driven into the more massive cluster is stronger on the side from which the less massive cluster approaches. However, 
at 0.5 Gyr prior to maximum luminosity the right-hand shock has not yet penetrated the more massive cluster core, as it has by this 
point in the equal-mass collisions. 

Peak luminosity and temperature occur as the core of the less massive cluster passes through (in the head-on case) or swings past 
(in the offset cases) that of the more massive cluster As in the equal-mass runs, in the offset cases here an S-shaped region of hot gas 
is apparent in the region between the cluster cores. The heating of the gas is more severe in the region in front of the less massive 
cluster Also, the less massive cluster undergoes extreme distortions at this time. Gas velocities of up to 2000 km s"' are present 
behind it, and as it plows into the hotter, denser ICM of the larger cluster, it is flattened perpendicular to its direction of motion. The 
core of the larger cluster is less affected by the interaction, as expected. Shortly after first core passage, the offset case C5 resembles 
Abell 754, consistent with the results of Roettiger et al. (1998). The offset cases also resemble Abell 2256 just at core passage if one 
uses the ASCA or BeppoSAX temperature maps (Markevitch 1996; Molendi, de Grandi, & Fusco-Femiano 2000) and assumes the 
merger axis to be perpendicular to the plane of the sky. 

Although the smaller cluster undergoes severe distortion during the core passage, the merger's effect on the total luminosity and 
average emission-weighted temperature of the system is noticeably less than in the 1 : 1 collisions. The luminosity increases by a factor 
of 2-4, and the temperature by a factor of 1 .4-2.2, for about 1/2 sound-crossing time during the initial core passage. While maximum 
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Fig. 6. — Surface density maps for dark matter (colormap) and gas (contours) in the 1:3 collision runs (C4— C6), projected along the z-axis. Scaling of the plot is 
the same as in Figure 3. 
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Fig. 7. — 2-10 keV emission-weighted temperature Tx (colormap) and X-ray surface briglitness Sx (contours) in tlie 1:3 collision runs (C4-C6), as seen along 
the z-axis. Scaling of the plot is the same as in Figure 4. 




-20246 
Time {t-t^J/t^^ 

Fig. 8. — 2—10 keV luminosity Lx and average emission-weighted temperature Tx as functions of time 1 in the 1:3 collision runs (C4-C6). Scahng of the plot is 
the same as in Figure 5. 
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Fig. 9. — Gas specific entropy and velocity fields in the collision plane (xv") of each run ait = 14.5 Gyr, near the end of the simulations. Entropy (colormap) is 
measured relative to the initial central value of the more massive cluster in each run in units of 1. 5k. The color key to the right of the figure indicates the range of 
values for this variable. Gas velocity is indicated by arrows; the fiducial arrow at the lower left has length 560 km s"'. The length scale bar at the upper left has 
length 0.5/i"' Mpc. 

separation following core passage is reached at about the same time as in the 1 : 1 cases, the drop in luminosity and temperature is 
not as great. This period is again followed by a weaker secondary peak in luminosity and temperature as the cores fall back together 
Here the most offset case (run C6) spends more time at maximum separation, causing the spread in secondary peak times to be 50% 
greater for the 1:3 runs than in the 1:1 cases. The second core passage for this run takes place almost 7 Gyr after the first, and in 
all cases this second interaction takes place on a nearly radial trajectory due to the dissipation caused by the first interaction. (The 
last two panels of Figures 6 and 7 for this run do not show the less massive cluster because at these times it is at a separation which 
places it outside the plotted region.) Both luminosity and temperature settle down to final values of 50% and 1 10%, respectively, of 
their initial values following the second core passage. A few oscillations are again seen as the cores settle into a common remnant. 
Unlike the 1 : 1 cases, here there is little variation in the final luminosity as impact parameter is varied. 



4. DISCUSSION 
4. 1 . Structure of Merger Remnants 

Each merger simulation was followed for approximately 15 Gyr, or about six sound-crossing times after the first core passage. 
Although this time is greater than the age of the Universe for il = 1 and h = 0.6, it is nevertheless insufficient time for the remnants of 
either the 1:1 or the 1:3 mergers to come to equilibrium. This point is illustrated by Figure 9, which shows the gas specific entropy 
and velocity fields in the xy-plane of each run at f = 15 Gyr Although in each case the core of the remnant has lower entropy than its 
surroundings, the gas at the margins (r > 5rs) continues to be disturbed by small-scale convective motions as matter accretes onto the 
remnant behind the shocks that have been driven out of the core region. The offset mergers all show evidence of bulk rotation, with 
velocities of ^ 200-300 km s"', but in even the most offset cases this provides at most a few percent of the support needed to keep 
the remnant in hydrostatic equilibrium. The 1 : 1 cases are closest to equilibrium, having settled into roughly spherical distributions 
within r ^ Rof the center of mass. The 1 :3 cases are relatively further from equilibrium, with distinct streams of material falling in 
from the right (in the head-on case) and along the top edge (in the offset cases). These streams continue to feed high-entropy gas 
from the smaller cluster into the remnant as the simulations come to an end. 

Departures from equilibrium can also be described in terms of total energetics. In describing the results of gas-only merger 
simulations, Ricker (1998) defined a 'virial disequilibrium' parameter as 



y 



2(T + U) 



W 



+ 1 



(36) 



where T, U, and W are the total kinetic, internal, and potential energies. When the system is in virial equilibrium, V = 0, and when it 
is not, V is equal to the second time derivative of the moment of inertia, divided by W. Thus a useful equilibrium criterion is to require 
y < e < < 1 for an extended period (say ~ fjc)- In the gas-only mergers, the clusters required between 5 .5 and 6 sound-crossing times 
to reach the level e r^ 0.02, after which V was roughly constant. In Figure 10 we plot V versus time for the six gas-plus-dark matter 
mergers considered in this paper Note that, as in the earlier work, V passes through a maximum during the first core passage, when 
the system is furthest from equilibrium. V then declines steeply but oscillates for several sound-crossing times as the dark matter 
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Fig. 10. — Virial disequilibrium parameter in the six collision runs. Time is offset from the time of peak luminosity t„, 



, and scaled to the sound-crossing time fsc 



continues to pump the gas. Finally the system begins to tend steadily toward V = around 4fsc after first core passage. However, by 
the end of the simulations only runs C4 and C5 have satisfied the equilibrium criterion just described. 

We plot angle-averaged profiles of various quantities from the merger runs atf = 15 Gyr in Figures 11 and 12. Profiles are generated 
from the raw simulation data as for the single-cluster test runs described in § 3. 1 . Although the merger remnant is not in hydrostatic 
equilibrium, the average profiles still serve to characterize the effect of the collision on the structure of the clusters. As for the single- 
cluster runs, the profiles here are scaled to the initial values of p,, r,, and so forth, with the exception of the enclosed gas fraction, the 
gas radial Mach number, and the pressure support ratio. For the 1:3 mergers the scaling parameters used are those pertaining to the 
more massive cluster (B). 

The dark matter density profile (Figures 11a and 12a) shows very little change in shape as a result of the collisions. Of all of 
the properties of the merger remnant, its shape is also the least affected by variations in mass ratio and impact parameter The 
normalization (^ p, since dark matter dominates the potential) increases in each case because the merger remnant is more massive 
than the individual colliding clusters, but the scale radius changes very little. The asymptotic slope of the profile does not increase as 
in the single-cluster test runs. 

Figures lib and 12b give the profile of the dark matter velocity dispersion at the end of the simulations. At the end of the 1:1 
collisions, the central velocity dispersions are all increased as a result of the deeper potential well of the merger remnant. A similar 
but smaller effect is seen in the 1:3 mergers. At large radii, the velocity dispersions are all small as a result of the evaporation of 
faster moving particles and the outflow-only outer boundary conditions in the simulations. A similar effect was seen in the single 
cluster simulations (Figure 2b). 

Although the dark matter density profiles are not strongly affected by the mergers, the gas density profiles (Figures lie and 12c) 
change substantially as a result of the mergers. The central gas density in each case actually shows a substantial decrease from its 
initial value, dropping to 30-70% of its initial value for the 1:1 runs and to 20-30% of its initial value for the 1:3 runs. Using the 
half-density radius as an estimate for the core radius, we find that r^ increases by 20-25%. The final value of pgo shows a significant 
dependence on impact parameter For the equal-mass collisions, the final gas density is about twice as high for head-on mergers 
as for the most offset mergers. For the 1:3 collisions, the variation is only 50%, but the sense is reversed, with head-on collisions 
resulting in a lower central gas density. In all cases we see a steepening of the asymptotic slope from r~^ to r"^, as in the single-cluster 
runs, because of the initial cutoff at r = R combined with the outflow boundary conditions. 

The flattening of the gas density profile occurs because the shocks generated by the merger add entropy to the system and increase 
the buoyancy of the ICM in the relatively constant dark matter potential. The resulting profile for the enclosed gas fraction (Figures 
lid and 12d) typically begins below the initial value of ^^ 1% for radii comparable to the zone spacing, then increases roughly as 
^0.8-1.0 ^j^j almost meets the initial profile at r '--^ R/2 before flattening out to an asymptotic value slightly below the initial value. 
The change in curvature at small radii is also seen in the single-cluster runs and is a numerical effect due to the influence of particle- 
mesh force smoothing on the cuspy dark matter profile. Likewise, the drop from the initial profile between r ^ R/2 and r ~ /J is a 
numerical effect due to the expansion of material beyond the imposed initial density cutoff. The drop in the asymptotic value of the 
gas fraction occurs because more gas than dark matter is lost through the outflow boundary - ^ 30% versus ~ 16% - as shocks tend 
to eject collisional material while leaving the collisionless dark matter alone. The addition of an external confining pressure would 
most likely allow the asymptotic gas density slope to be flatter but leave the inner core unchanged. 

This is apparent from the gas temperature profiles (Figures lie and 12e), which show an inner region whose behavior depends 
on impact parameter and mass ratio and an outer region that is roughly independent of impact parameter and corresponds to the 
single-cluster result. The boundary between these two regions occurs at r ^ 3rc in each of the runs and is ^ 20% hotter relative to 
Tq in the 1:1 runs than in the 1:3 runs. In all cases the final temperature profile declines more slowly at the very center than does the 
initial profile. Offset mergers display opposite trends with increasing b for the two mass ratios: for 1:1 mergers, increasing impact 
parameter increases the final temperature over the head-on case, while for 1:3 mergers it decreases the temperature, reproducing the 
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Fig. 11. — Average profiles of dark matter and gas quantities at ? = 15 Gyr in the 1:1 collision runs (C1-C3). In all cases except (d), (g), and (h), quantities are 
scaled to their values at the beginning of the simulation, (a) Dark matter density, (b) Dark matter velocity dispersion, (c) Gas density, (d) Enclosed gas fraction, (e) 
Gas temperature, (t) Gas specific entropy, (g) Gas radial Mach number (h) Ratio of pressure force to gravitational force on the gas. 



initial profile in the most offset 1:3 case. This behavior occurs because of the increase in core entropy associated with mixing of 
shocked gas from the outer ICM into the core material, as discussed in § 4.3. The specific entropy profile (Figures 11 f and 12f) also 
demonstrates a parameter-independent outer region and parameter-dependent core, with the central entropy increase being ^ \.5k 
over the initial central value in the hottest runs. Note that the magnitude of this increase is roughly the same for the two mass ratios. 
Although the gas has not reached virial equilibrium in the simulations by f = 15 Gyr, it is relatively quiescent. Figures 11 g and 12g 
show the radial Mach number profiles for each run. Within r = lOrc the radial Mach number is less than 0.06 in magnitude, indicating 
that radial motions near the center are fairly subsonic at this point. Outside r = lOr^, the magnitude of the Mach number increases 
rapidly, reflecting the continued accretion of gas that was previously ejected from the center by the merger shocks. The gas within 
the initial cutoff radius is also quite close to hydrostatic equilibrium in angle average. The angle-averaged ratio of the pressure force 
to gravitational force per unit mass is shown in Figures 1 Ih and 12h. Except for the innermost few zones where the small number of 
zones in the radial bins and truncation error in the gradient calculation dominate this ratio, it is within 15% of unity for r < lOr^. 
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4.2. Luminosity and Temperature Evolution 

In the simplest self-similar picture of their assembly, clusters undergo spherical collapse and virialization and then radiate via 
thermal bremsstrahlung (Kaiser 1986). This model fails to predict the observed X-ray luminosity-temperature relationship: the 
predicted relation is Lx oc r|, whereas the observed relation is closer to Lx oc T^ (Mushotzky 1984). Cooling flows have larger 
luminosities and lower average temperatures than average and are known to steepen this relation and increase its dispersion (Fabian 
et al. 1994). However, removing them from cluster samples in various ways does not eliminate the discrepancy (Markevitch 1998; 
Allen & Fabian 1999; Arnaud & Evrard 1999). Various solutions to this problem have been proposed. Most recently these have 
focused on raising the entropy of the ICM via nongravitational heating, though Bryan (2000) has recently suggested that galaxy 
formation may do the same thing by removing low entropy gas. Preheating models increase the entropy of the gas before it falls 
into cluster potentials (Kaiser 1991; Evrard & Henry 1991). These models require substantial amounts of heating: ^ 1 keV/baryon 
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appears to be needed (Lloyd-Davies, Ponman, & Cannon 2000). 

Note that the luminosity-temperature relation used to scale the initial conditions for our simulations (eq. [24]) is not retained as 
the mergers progress. Comparing T^^ to Lx as time passes in run CI, we find that immediately following the initial brightness peak, 
the ratio of these two quantities jumps to twice its initial value, fluctuates here for a bit, increases briefly to almost six times its initial 
value during the second core passage, then settles down to about three times its initial value, where it remains roughly constant for 
the remainder of the calculation. The behavior of the other runs is similar Thus the trajectory of the system in the Lx - Tx plane is 
roughly parallel to the observed statistical relation both before and after the merger, but during the merger the system undergoes a 
shift to a higher temperature normalization. 

The normalization of the Lx - Tx relation is complicated by the presence of two clusters within the computational volume. In this 
context we note that luminosity is an 'extensive' quantity, which must be divided by the number of clusters one is counting, while 
temperature is an 'intensive' quantity, which does not scale with the number of clusters. For equal-mass mergers, the cluster cores 
merge to form a single core (to within the spatial resolution of our grid) approximately 5-7 Gyr after the initial core passage. Long 
before this point the clusters would be considered a single cluster with 'bimodal' structure if seen along the z-axis, because the twin 
brightness peaks lie within an Abell radius of each other and are gravitationally bound. Moreover, from most lines of sight the cluster 
would appear to be unimodal by at most 3 Gyr after initial core passage. Indeed, the head-on collision seen along the collision axis 
appears to be unimodal throughout the calculation. Comparing the initial and final states of the system is less problematic. Initially 
we have two clusters following the imposed Lx - Tx relation; at the end, we have one cluster with temperature 1 .4 times the initial 
value and luminosity between 80% and 160% of the initial value of the individual clusters. Hence the final ratio of T^'' to Lx is in 
the range L5 (run CI) to 3 (run C3). 

Alternatively, one may say that the final temperature is ^ L2- L5x the value expected given the final luminosity and the initial 
Lx-Tx relation. Let us define 



f^final /^initial 



'i/r^™ 



/'r final / HnitiaA V^^i 



(37) 



In the highest-resolution single-cluster run (S7) we see fir ^ 1.04 due to numerical effects. Nearly all of the X-ray emission arises 
within a few scale radii of the center of the cluster, so the increase in run S7 (and presumably a comparable portion of the increase in 
the merger runs) is probably due to spatial resolution. The lack of an external confining pressure (such as might be expected due to 
continuously infalling material from the low-density intergalactic medium) allows material from the center of the cluster to expand 
more than it might otherwise, causing unwanted adiabatic cooling. Although this changes fir, adiabatic changes yield Lx oc 7^, so 
cooling due to the outflow boundaries should decrease this ratio, not increase it: 
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(38) 



Shocks might have a much larger effect on fu because they produce less of a density jump than a temperature jump. For shocks of 
arbitrary strength, the inverse of the shock compression is given by (Markevitch, Sarazin, & Vikhlinin 1999): 
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where px {pi) and Tx {Ti) are the preshock (postshock) density and temperature, respectively, and 7 is the adiabatic index. For 
7 = 5/3, this gives 
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(40) 



The average temperature increases by about 40% during the equal-mass merger runs, so the value of fu lies between L2 and L5. 
In the 1:3 cases, the less massive cluster contributes only about 10% as much luminosity as the more massive cluster at the start of 
the calculation, so in this case fu ~ 14, consistent with the range of values in the 1:1 collisions. On the other hand, equations (38) 
and (40) give /l7-| adiabatic ~ 1-08 and /l7-| shock ~ 1.10. Thus, the temperature increase is larger (or the X-ray luminosity increase is 
smaller) than expected from either shocks or adiabatic compression. This probably indicates that the luminosity increase is limited 
by the finite numerical resolution in the cores of the clusters and would be larger if the calculations were better resolved. 

We close this section with the following note. The observation of massive clusters at high redshift (e.g., MS 1054: Donahue et al. 
1998) is often taken as strong evidence for low J7o because of the exponential sensitivity of the cluster mass function to i^o (e-g-, 
Oukbir & Blanchard 1992). However, in deriving such constraints one directly measures an X-ray temperature, then uses the mass- 
temperature relation to find the mass. Our results show that this procedure is dangerous: merger shocks produce short-lived, but 
significant, increases in the X-ray temperature of a cluster. For a few percent of its life, a merging cluster appears much hotter, and 
thus more massive, than it really is (Figure 13). Furthermore, because this period corresponds to a time of significantly enhanced 
luminosity, merging clusters are also easier to detect. We note, for example, that MS 1054 now shows evidence of substructure and a 
lower intrinsic temperature than first reported, suggesting that its apparent high temperature is due to an ongoing merger (Neumann 
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& Arnaud 2000). As another example, Gomez, Hughes, and Birkinshaw (2000) have recently suggested that Abell 665, one of the 
richest, hottest, and most X-ray luminous clusters at redshifts z < 0.2, is currently undergoing a merger and is at the epoch when the 
cluster cores collide. This could explain its very high luminosity and temperature. We are using our results with a merger tree code 
to estimate the effect of merger shocks on the significance of distant, massive cluster constraints on VIq and will address this question 
in a subsequent paper (Randall, Sarazin, & Ricker 2001, in preparation). 

4.3. Persistence of Gas Cores and Turbulence 

A cluster is defined by the presence of a distinct potential well within which dark matter, gas, and galaxies are confined. Mergers 
produce dramatic fluctuations in these potentials, leading eventually to the combination of the merging clusters into a single rem- 
nant. Some merging clusters in which the mergers appear to be well-advanced show a double core structure at their centers (e.g., 
Abell 2065; Markevitch et al. 1999), with two X-ray surface brightness peaks, indicating that the cores of the individual cluster 
potentials are still intact. If the original subclusters have cooling flows at their centers, the X-ray surface brightness of the surviving 
peaks in the subcluster potentials can be particularly large. 

Cluster cores lose their distinct identity in several ways. As with galaxies that pass too close to a cluster's center, tidal stripping can 
remove material from a shallower cluster potential and add it to a deeper one. In addition, violent relaxation (for collisionless matter) 
and ram-pressure stripping (for diffuse gas) convert the bulk kinetic energy of a cluster into dispersive or thermal energy, enabling 
material to become unbound from its own potential, while remaining bound within the combined potential of the material from both 
clusters. Dynamical friction offers another means for converting organized bulk velocities to dispersive ones, but the timescale on 
which it operates is typically large compared to the timescales associated with these other mechanisms. 

Our results shed some light on the mechanism by which ram pressure affects the survival of cluster cores in mergers. Rather than 
ablating gas from the edges of a core, ram pressure disrupts cores by pushing the gas as a unit away from the center of the potential 
well in which it originates, causing it to become convectively unstable. The resulting convective motions destroy the spherical 
symmetry of the core and mix the core gas with higher-entropy shocked gas in its surroundings. This also appears to be an important 
mechanism by which mergers initiate turbulent motions. 

Figure 14 shows an example of this process using the specific entropy field. This example begins just after the first core interaction 
of run C5 (the less massive cluster core has exited to the right at 4.00 Gyr). The diffusivity of the gas is low, so the core interaction 
pushes low-entropy gas out of the more massive core, where it finds itself out of convective equilibrium. Beginning around 5.00 Gyr 
this gas flows in a plume back toward the more massive potential center. A large eddy also forms in the wake of the less massive 
cluster; its remnants persist until 7.00 Gyr The plume overshoots the potential center and begins to decelerate and spread along 
equipotential lines. It reaches maximum overshoot around 7.00 Gyr, when it is met by high-entropy gas from the returning less 
massive core. The remaining gas bound to this core has higher entropy and, stripped from its potential well, is forced to flow around 
the plume. As the dark matter associated with the less massive cluster passes through the more massive potential well, it drives a 
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Fig. 14. — A merger-driven convective plume. The plots show specific entropy, measured relative to the initial central value of the more massive cluster, in the 
collision plane of run C5. The colormap is chosen to accentuate low-entropy regions, so the shocked region in the first four frames, where the colormap is saturated, 
appears white. Contours show the gravitational potential, while arrows show the gas velocity field. The fiducial bar at lower left has a length of 250/i"' kpc. 



new shock through the low-entropy gas from the more massive cluster that did not flow into the plume; this shock exits to the lower 
left after 8.50 Gyr The remnants of the plume mix with the newly reshocked gas to form the core of the merger remnant, with higher 
entropy than was present in the original clusters. 

The turbulent velocity field created by a merger lasts for several sound-crossing times. In Figure 15 we plot the rms gas velocity 
profile in each of the collision runs at f w 15 Gyr as a fraction of the average sound speed and the circular (Keplerian) velocity. 
Throughout the merger remnant, turbulent velocities are subsonic: for r < 2Qrc they are typically 10-20% of the sound speed in the 
1:1 cases, while the rms Mach number is about twice as large in the 1:3 cases. Outside this radius the rms velocity rises to between 
60% and 100% of the sound speed, with the increase much steeper in the 1:1 cases. Throughout most of the cluster, turbulence 
supplies '--^ 5 - 10% of the hydrostatic support. In certain cases this contribution appears to be larger within the innermost r^. or so 
of the core, though it is likely that here the rms velocity is overestimated due to the relatively small number of zones contributing to 
radial velocity averages. 

Our results support the finding of others (e.g., Roettiger, Burns, & Loken 1996; Norman & Bryan 1998) that, following a merger, 
the intracluster medium becomes turbulent, and that turbulent pressure and bulk rotation provide less than 20% of the support 
needed for hydrostatic equilibrium. The picture which emerges from our calculations is that this turbulence is initiated by convective 
instability and, to a lesser extent, by baroclinic vorticity generation at oblique shocks. It is sustained for several Gyr by pumping 
due to the oscillation of the dark matter, which relaxes more slowly than the gas. While Roettiger et al. (1996) also argue for this 
mechanism using isolated merger simulations, Norman and Bryan (1998) find higher rms velocities in a cluster simulation involving 
multiple mergers and accretion, and they argue that minor mergers are the primary pumping mechanism. Our results suggest that 
individual major mergers and cumulative minor mergers induce comparable amounts of turbulence, since the turbulent pressure in 
our simulations is ~ 1 /2 - 1 /3 that of Norman & Bryan's result even several Gyr after the collision. 

To date the grid-based merger simulations used to argue for turbulence in clusters have reached spatial resolutions on the order of 
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10 kpc. In fully developed turbulence, the range of eddy length scales is of order 

r-(R^)'^\ (41) 

where £q is the integral scale and £4 is the dissipation scale (Richardson 1922; Kolmogoroff 1941). Thus, directly simulating three- 
dimensional turbulent flow with a hydrodynamic code requires at least of order {Refl"' zones. For many astrophysical applications, 
this is impossible - for example, in stellar convective layers Re ^ 10^"*. However, as we have already argued, the Reynolds number 
of the intracluster medium during mergers is of order 10^, and turbulence, if it is present at all, operates over a much smaller range of 
length scales. We therefore have some hope of performing direct numerical simulations of merger-driven turbulence within the very 
near future. By comparison, the simulations described in this paper achieve a spatial dynamic range of ^ 200. 

Several approximate methods for numerically treating turbulence exist (see Canuto 1994 for a partial review). Among these meth- 
ods, large-eddy simulation (LES) has achieved some success in reproducing characteristics of turbulence in laboratory experiments 
(Canuto 1997). LES attempts to simulate directly only the largest-scale eddies in a turbulent flow, since these depend on the bound- 
ary conditions and geometry of the flow at the integral scale, and it uses a subgrid model to handle the (presumably universal and 
geometry-independent) turbulent dissipation on unresolved scales. Most grid-based astrophysical simulation algorithms are, on some 
level, LES, but they do not use subgrid turbulence models, instead relying on numerical dissipation at the zone scale. This numerical 
dissipation arises due to the truncation error of each individual method and bears no relation to the real dissipation mechanism. For 
most methods this practice should yield an incorrect energy dissipation spectrum, since it links large- and small-scale eddies in an 
unphysical way (Canuto 1997). However, it has been argued by Sytine et al. (2000) that the piecewise-parabolic method, which 
we employ in this paper, accurately reproduces the energy and enstrophy power spectra produced by solving the full Navier-Stokes 
equations in direct simulations of compressible turbulence, because it is much less dissipative than other techniques and because its 
dissipation is confined to narrow regions near flow discontinuities. Estimates of the dissipation in PPM for some simple flow fields 
(Porter & Woodward 1994) suggest that convective features spanned by 32 or more zones are well-converged. This is consistent with 
recent results of Calder et al. (2001) for single-mode Rayleigh-Taylor instabilities. Thus, for example, the convective plume shown in 
Figure 14, which measures approximately 32x 10 zones at maximum extent, is probably just sufficiently resolved to yield the correct 
growth rate. Further high-resolution calculations using adaptive mesh refinement will be useful in verifying this result (Ricker et al. 
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2001). 



5. SUMMARY 



We have presented results from a controlled parameter study of off-axis mergers between clusters of galaxies. We have extended 
previous work by using the NFW density profile and nonisothermal temperature profiles, systematically controlling merger param- 
eters and justifying their scaling, carefully studying numerical effects and convergence requirements, and examining the survival of 
pressure peaks and the mechanism for the onset of turbulence in mergers. We have also presented results on the variation of relaxation 
times, merger remnant properties, and luminosity and temperature peaks with merger parameters. The calculations presented here 
do not include radiative cooling or magnetic fields. We anticipate adding these physical processes in future calculations (particularly 
cooling, which appears to play a very important dynamical role in many observed clusters). 

The morphological changes, relative velocities, and temperature jumps we observe agree well with previous studies of collisions 
between clusters modeled using the King profile (e.g., Roettiger et al. 1998). In particular, velocities of ^ 2000 km s~' are to be 
expected, even in quite offset mergers, for cosmologically reasonable initial conditions. The peak temperature occurs for a period 
< /sc/2 at first core passage and is a factor ^2-3 greater than the initial temperature. We observe a larger jump in X-ray luminosity 
(~ 4- lOx) than previous studies including dark matter, and we argue that this increase is most likely a lower limit due to our spatial 
resolution. Second core passage, with corresponding secondary luminosity and temperature jumps, typically occurs 1.5-2.5/sc after 
the first. In the most offset 1:3 case, this secondary peak is delayed by an additional sound-crossing time. We emphasize that 
luminosity and temperature jumps due to mergers may have an important bearing on constraints on fl derived from the observation 
of hot clusters at high redshift. We will address this issue further in a subsequent paper. 

Shocks play an important dissipative role in mergers, but they are relatively weak in the highest-density regions. As a result 
they do not directly raise the entropy of the cluster cores. Instead, shocks create entropy in the outer parts of the clusters, and this 
high-entropy gas is mixed with the core gas during later stages of the merger. Mixing is initiated by ram pressure: the core gas is 
displaced from its potential center and becomes convectively unstable. The resulting convective plumes initiate large-scale turbulent 
motions with eddy sizes up to several 100 kpc. This turbulence is pumped by oscillations in the gravitational potential, which in turn 
are driven by the more slowly relaxing collisionless dark matter. Even after nearly a Hubble time these motions persist as subsonic 
turbulence in the cluster cores, providing 5- 10% of the support against gravity. The dark matter oscillations are also reflected in 
the extremely long time following a merger required for the remnant to come to virial equilibrium. Because of the increase in core 
entropy, if a constant-density gas core is initially present, it is not destroyed by a merger, even if the dark matter density profile has a 
central cusp. 

This research has been supported in part by NASA Astrophysical Theory Program grant NAG 5-3057 and Chandra Award Num- 
bers GO0-1019X, GO0-1141X, GO0-1158X, and GO0-1173X at the University of Virginia, and by the ASCI Flash Center at the 
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burgh Supercomputing Center and the San Diego Supercomputer Center. P. M. R. would like to acknowledge useful conversations 
with D. Lamb, M. C. Miller, J. Stone, and E. Ostriker. 
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